##### PACKAGES #####
library(multiwayvcov)
library(lmtest)
library(dplyr)
library(ggplot2)

##### Data: isq #####

#Note: All analyses include fixed effects for parish and clustered standard errors for enumeration districts. 
#This was done to account for our sampling design, which blocked by parish and clustered by ED


### Creating the aspiring migrant subsample ###
potmig <- subset(isq, aspire.abroad.pre>4)
potmig$ed.parish <- paste(potmig$ed.1, potmig$parish, sep=".")



### Table 1 ###

### Full sample models
### Model 1: enter (illegal migration)
enter.ate <- lm(enterlist.pre ~ enterlist.pre.dum + factor(parish), data=MD.LIST)
summary(enter.ate)
enter.ate.clust<-cluster.vcov(enter.ate, MD.LIST$ed.parish)
enter.ate.mcse <- coeftest(enter.ate, enter.ate.clust)
enter.ate.mcse

### Model 3: work (semi-legal migration)
work.ate <- lm(worklist.pre ~ worklist.pre.dum + factor(parish), data=MD.LIST)
summary(work.ate)
work.ate.clust<-cluster.vcov(work.ate, MD.LIST$ed.parish)
work.ate.mcse <- coeftest(work.ate, work.ate.clust)
work.ate.mcse


### Aspiring migrant models

# Because the treatment was randomized at the level of the full sample,
# we include our control variables (gender, age, education, income, networks abroad) 

### Model 2: enter (illegal migration)
enter.ate.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + gender + age + education + income.imp + freq.help.net+  factor(parish), data=potmig)
summary(enter.ate.potmig)
enter.ate.potmig.clust<-cluster.vcov(enter.ate.potmig, potmig$ed.parish)
enter.ate.potmig.mcse <- coeftest(enter.ate.potmig, enter.ate.potmig.clust)
enter.ate.potmig.mcse

# Model 4: work (semi-legal migration)
work.ate.potmig <- lm(worklist.pre ~ worklist.pre.dum + gender + age + education + income.imp + freq.help.net+  factor(parish), data=potmig)
summary(work.ate.potmig)
work.ate.potmig.clust<-cluster.vcov(work.ate.potmig, potmig$ed.parish)
work.ate.potmig.mcse <- coeftest(work.ate.potmig, work.ate.potmig.clust)
work.ate.potmig.mcse

### Figure 1 ###

### Family and enter (illegal migration)
## running the model
enter.family.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + norm.family.rev + enterlist.pre.dum*norm.family.rev
                        + gender + age + education + income.imp + freq.help.net + factor(parish), 
                        data= potmig)

summary(enter.family.potmig)

##generating clustered standard errors
enter.family.potmig.clust<-cluster.vcov(enter.family.potmig, potmig$ed.parish)
enter.family.potmig.mcse <- coeftest(enter.family.potmig, enter.family.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.family.potmig.b2 <- enter.family.potmig.mcse[2,1]
enter.family.potmig.b3 <- enter.family.potmig.mcse[3,1]
enter.family.potmig.b22 <- enter.family.potmig.mcse[22,1] # interaction term
# variances
enter.family.potmig.varb2 <- enter.family.potmig.clust[2,2]
enter.family.potmig.varb3 <- enter.family.potmig.clust[3,3]
enter.family.potmig.varb22 <- enter.family.potmig.clust[22,22]
# covariances
enter.family.potmig.covb2b22 <- enter.family.potmig.clust[2,22]
enter.family.potmig.covb3b22 <- enter.family.potmig.clust[3,22]

list(enter.family.potmig.b2, enter.family.potmig.b3, enter.family.potmig.b22, enter.family.potmig.varb2, 
     enter.family.potmig.varb3, enter.family.potmig.varb22, enter.family.potmig.covb2b22, enter.family.potmig.covb3b22)

##transforming dataframe to graph    
enter.family.potmig.df <- potmig %>% mutate(norm.family.rev = seq(1, 7, length.out=802),
                                          enter.family.potmig.conbx = enter.family.potmig.b2+enter.family.potmig.b22*norm.family.rev,
                                          enter.family.potmig.consx = sqrt(enter.family.potmig.varb2 + enter.family.potmig.varb22*(norm.family.rev^2)+
                                                                           2*enter.family.potmig.covb2b22*norm.family.rev),
                                          enter.family.potmig.ax = 1.96*enter.family.potmig.consx,
                                          enter.family.potmig.upperx = enter.family.potmig.conbx + enter.family.potmig.ax,
                                          enter.family.potmig.lowerx = enter.family.potmig.conbx - enter.family.potmig.ax)


### Family and work (semi-legal migration)
## running the model

work.family.potmig <- lm(worklist.pre ~ worklist.pre.dum + norm.family.rev + worklist.pre.dum*norm.family.rev
                       + gender + age + education + income.imp + freq.help.net + factor(parish), 
                       data= potmig)

summary(work.family.potmig)

##generating clustered standard errors
work.family.potmig.clust<-cluster.vcov(work.family.potmig, potmig$ed.parish)
work.family.potmig.mcse <- coeftest(work.family.potmig, work.family.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.family.potmig.b2 <- work.family.potmig.mcse[2,1]
work.family.potmig.b3 <- work.family.potmig.mcse[3,1]
work.family.potmig.b22 <- work.family.potmig.mcse[22,1] # interaction term
# variances
work.family.potmig.varb2 <- work.family.potmig.clust[2,2]
work.family.potmig.varb3 <- work.family.potmig.clust[3,3]
work.family.potmig.varb22 <- work.family.potmig.clust[22,22]
# covariances
work.family.potmig.covb2b22 <- work.family.potmig.clust[2,22]
work.family.potmig.covb3b22 <- work.family.potmig.clust[3,22]

list(work.family.potmig.b2, work.family.potmig.b3, work.family.potmig.b22, work.family.potmig.varb2, 
     work.family.potmig.varb3, work.family.potmig.varb22, work.family.potmig.covb2b22, work.family.potmig.covb3b22)

## transforming dataframe to graph    
work.family.potmig.df <- potmig %>% mutate(norm.family.rev = seq(1, 7, length.out=802),
                                         work.family.potmig.conbx = work.family.potmig.b2+work.family.potmig.b22*norm.family.rev,
                                         work.family.potmig.consx = sqrt(work.family.potmig.varb2 + work.family.potmig.varb22*(norm.family.rev^2)+
                                                                         2*work.family.potmig.covb2b22*norm.family.rev),
                                         work.family.potmig.ax = 1.96*work.family.potmig.consx,
                                         work.family.potmig.upperx = work.family.potmig.conbx + work.family.potmig.ax,
                                         work.family.potmig.lowerx = work.family.potmig.conbx - work.family.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

## The rug
family.potmig.temp <- potmig$norm.family.rev
a<-runif(802,min=0,max=0.6)
a<-a-0.2
family.potmig.jitter<-family.potmig.temp+a

## The plot
family.potmig.df <- cbind(work.family.potmig.df, enter.family.potmig.df)
colnames(family.potmig.df) <- make.unique(names(family.potmig.df))

family.potmig.plot.combined <- ggplot(data = family.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(1,7)) +
  geom_ribbon(aes(x=norm.family.rev, ymin=work.family.potmig.lowerx, ymax = work.family.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=norm.family.rev, y=work.family.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=norm.family.rev, ymin=enter.family.potmig.lowerx, ymax = enter.family.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=norm.family.rev, y=enter.family.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=family.potmig.jitter),sides="b") +
  xlab("OK to Break Law to Help Family") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+ 
  theme(legend.title = element_blank())

family.potmig.plot.combined


### Figure 2 ###

### Unfairness and enter (illegal migration)
## running the model

enter.fair.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + fair.jam+ enterlist.pre.dum*fair.jam
                        + gender + age + education + income.imp + freq.help.net + factor(parish), 
                        data= potmig)

summary(enter.fair.potmig)

## generating clustered standard errors
enter.fair.potmig.clust<-cluster.vcov(enter.fair.potmig, potmig$ed.parish)
enter.fair.potmig.mcse <- coeftest(enter.fair.potmig, enter.fair.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.fair.potmig.b2 <- enter.fair.potmig.mcse[2,1]
enter.fair.potmig.b3 <- enter.fair.potmig.mcse[3,1]
enter.fair.potmig.b22 <- enter.fair.potmig.mcse[22,1] # interaction term
# variances
enter.fair.potmig.varb2 <- enter.fair.potmig.clust[2,2]
enter.fair.potmig.varb3 <- enter.fair.potmig.clust[3,3]
enter.fair.potmig.varb22 <- enter.fair.potmig.clust[22,22]
# covariances
enter.fair.potmig.covb2b22 <- enter.fair.potmig.clust[2,22]
enter.fair.potmig.covb3b22 <- enter.fair.potmig.clust[3,22]

list(enter.fair.potmig.b2, enter.fair.potmig.b3, enter.fair.potmig.b22, enter.fair.potmig.varb2, 
     enter.fair.potmig.varb3, enter.fair.potmig.varb22, enter.fair.potmig.covb2b22, enter.fair.potmig.covb3b22)

##transforming dataframe to graph    
enter.fair.potmig.df <- potmig %>% mutate(fair.jam = seq(1, 7, length.out=802),
                                          enter.fair.potmig.conbx = enter.fair.potmig.b2+enter.fair.potmig.b22*fair.jam,
                                          enter.fair.potmig.consx = sqrt(enter.fair.potmig.varb2 + enter.fair.potmig.varb22*(fair.jam ^2)+
                                                                           2*enter.fair.potmig.covb2b22*fair.jam),
                                          enter.fair.potmig.ax = 1.96*enter.fair.potmig.consx,
                                          enter.fair.potmig.upperx = enter.fair.potmig.conbx + enter.fair.potmig.ax,
                                          enter.fair.potmig.lowerx = enter.fair.potmig.conbx - enter.fair.potmig.ax)


### Unfairness and work (semi-legal migration)
## running the model

work.fair.potmig <- lm(worklist.pre ~ worklist.pre.dum + fair.jam + worklist.pre.dum*fair.jam
                       + gender + age + education + income.imp + freq.help.net + factor(parish), 
                       data= potmig)

summary(work.fair.potmig)

## generating clustered standard errors
work.fair.potmig.clust<-cluster.vcov(work.fair.potmig, potmig$ed.parish)
work.fair.potmig.mcse <- coeftest(work.fair.potmig, work.fair.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.fair.potmig.b2 <- work.fair.potmig.mcse[2,1]
work.fair.potmig.b3 <- work.fair.potmig.mcse[3,1]
work.fair.potmig.b22 <- work.fair.potmig.mcse[22,1] # interaction term
# variances
work.fair.potmig.varb2 <- work.fair.potmig.clust[2,2]
work.fair.potmig.varb3 <- work.fair.potmig.clust[3,3]
work.fair.potmig.varb22 <- work.fair.potmig.clust[22,22]
# covariances
work.fair.potmig.covb2b22 <- work.fair.potmig.clust[2,22]
work.fair.potmig.covb3b22 <- work.fair.potmig.clust[3,22]

list(work.fair.potmig.b2, work.fair.potmig.b3, work.fair.potmig.b22, work.fair.potmig.varb2, 
     work.fair.potmig.varb3, work.fair.potmig.varb22, work.fair.potmig.covb2b22, work.fair.potmig.covb3b22)

## transforming dataframe to graph    
work.fair.potmig.df <- potmig %>% mutate(fair.jam = seq(1, 7, length.out=802),
                                         work.fair.potmig.conbx = work.fair.potmig.b2+work.fair.potmig.b22*fair.jam ,
                                         work.fair.potmig.consx = sqrt(work.fair.potmig.varb2 + work.fair.potmig.varb22*(fair.jam ^2)+
                                                                         2*work.fair.potmig.covb2b22*fair.jam ),
                                         work.fair.potmig.ax = 1.96*work.fair.potmig.consx,
                                         work.fair.potmig.upperx = work.fair.potmig.conbx + work.fair.potmig.ax,
                                         work.fair.potmig.lowerx = work.fair.potmig.conbx - work.fair.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

# The rug
fair.potmig.temp <- potmig$fair.jam
a<-runif(802,min=0,max=0.6)
a<-a-0.2
fair.potmig.jitter<-fair.potmig.temp+a

# The plot
fair.potmig.df <- cbind(work.fair.potmig.df, enter.fair.potmig.df)
colnames(fair.potmig.df) <- make.unique(names(fair.potmig.df))

fair.potmig.plot.combined <- ggplot(data = fair.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(1,7)) +
  geom_ribbon(aes(x=fair.jam, ymin=work.fair.potmig.lowerx, ymax = work.fair.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=fair.jam, y=work.fair.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=fair.jam, ymin=enter.fair.potmig.lowerx, ymax = enter.fair.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=fair.jam, y=enter.fair.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=fair.potmig.jitter),sides="b") +
  xlab("Immigration Authorities are Unfair") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+ 
  theme(legend.title = element_blank())

fair.potmig.plot.combined

### Figure 3 ###

### Perceived risk and enter (illegal migration)
## running the model

enter.risk.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + risk.enter.dest + enterlist.pre.dum*risk.enter.dest 
                        + gender + age + education + income.imp + freq.help.net + factor(parish), 
                        data= potmig)

summary(enter.risk.potmig)

## generating clustered standard errors
enter.risk.potmig.clust<-cluster.vcov(enter.risk.potmig, potmig$ed.parish)
enter.risk.potmig.mcse <- coeftest(enter.risk.potmig, enter.risk.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.risk.potmig.b2 <- enter.risk.potmig.mcse[2,1]
enter.risk.potmig.b3 <- enter.risk.potmig.mcse[3,1]
enter.risk.potmig.b22 <- enter.risk.potmig.mcse[22,1] # interaction term
# variances
enter.risk.potmig.varb2 <- enter.risk.potmig.clust[2,2]
enter.risk.potmig.varb3 <- enter.risk.potmig.clust[3,3]
enter.risk.potmig.varb22 <- enter.risk.potmig.clust[22,22]
# covariances
enter.risk.potmig.covb2b22 <- enter.risk.potmig.clust[2,22]
enter.risk.potmig.covb3b22 <- enter.risk.potmig.clust[3,22]

list(enter.risk.potmig.b2, enter.risk.potmig.b3, enter.risk.potmig.b22, enter.risk.potmig.varb2, 
     enter.risk.potmig.varb3, enter.risk.potmig.varb22, enter.risk.potmig.covb2b22, enter.risk.potmig.covb3b22)

## transforming dataframe to graph    
enter.risk.potmig.df <- potmig %>% mutate(risk.enter.dest = seq(0, 10, length.out=802),
                                          enter.risk.potmig.conbx = enter.risk.potmig.b2+enter.risk.potmig.b22*risk.enter.dest,
                                          enter.risk.potmig.consx = sqrt(enter.risk.potmig.varb2 + enter.risk.potmig.varb22*(risk.enter.dest^2)+
                                                                           2*enter.risk.potmig.covb2b22*risk.enter.dest),
                                          enter.risk.potmig.ax = 1.96*enter.risk.potmig.consx,
                                          enter.risk.potmig.upperx = enter.risk.potmig.conbx + enter.risk.potmig.ax,
                                          enter.risk.potmig.lowerx = enter.risk.potmig.conbx - enter.risk.potmig.ax)


### Perceived risk and work (semi-legal migration)
## running the model

work.risk.potmig <- lm(worklist.pre ~ worklist.pre.dum + risk.work.dest + worklist.pre.dum*risk.work.dest 
                       + gender + age + education + income.imp + freq.help.net + factor(parish), 
                       data= potmig)

summary(work.risk.potmig)

## generating clustered standard errors
work.risk.potmig.clust<-cluster.vcov(work.risk.potmig, potmig$ed.parish)
work.risk.potmig.mcse <- coeftest(work.risk.potmig, work.risk.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.risk.potmig.b2 <- work.risk.potmig.mcse[2,1]
work.risk.potmig.b3 <- work.risk.potmig.mcse[3,1]
work.risk.potmig.b22 <- work.risk.potmig.mcse[22,1] # interaction term
# variances
work.risk.potmig.varb2 <- work.risk.potmig.clust[2,2]
work.risk.potmig.varb3 <- work.risk.potmig.clust[3,3]
work.risk.potmig.varb22 <- work.risk.potmig.clust[22,22]
# covariances
work.risk.potmig.covb2b22 <- work.risk.potmig.clust[2,22]
work.risk.potmig.covb3b22 <- work.risk.potmig.clust[3,22]

list(work.risk.potmig.b2, work.risk.potmig.b3, work.risk.potmig.b22, work.risk.potmig.varb2, 
     work.risk.potmig.varb3, work.risk.potmig.varb22, work.risk.potmig.covb2b22, work.risk.potmig.covb3b22)

## transforming dataframe to graph    
work.risk.potmig.df <- potmig %>% mutate(risk.work.dest = seq(0, 10, length.out=802),
                                         work.risk.potmig.conbx = work.risk.potmig.b2+work.risk.potmig.b22*risk.work.dest,
                                         work.risk.potmig.consx = sqrt(work.risk.potmig.varb2 + work.risk.potmig.varb22*(risk.work.dest^2)+
                                                                         2*work.risk.potmig.covb2b22*risk.work.dest),
                                         work.risk.potmig.ax = 1.96*work.risk.potmig.consx,
                                         work.risk.potmig.upperx = work.risk.potmig.conbx + work.risk.potmig.ax,
                                         work.risk.potmig.lowerx = work.risk.potmig.conbx - work.risk.potmig.ax)

### interaction plot with confidence intervals based on clustered SEs

# The rug
enter.risk.potmig.temp <- potmig$risk.enter.dest
work.risk.potmig.temp <- potmig$risk.work.dest
a<-runif(802,min=0,max=0.6)
a<-a-0.2
enter.risk.potmig.jitter<-enter.risk.potmig.temp+a
work.risk.potmig.jitter<-work.risk.potmig.temp+a

# The plot
risk.potmig.df <- cbind(work.risk.potmig.df, enter.risk.potmig.df)
colnames(risk.potmig.df) <- make.unique(names(risk.potmig.df))

risk.potmig.plot.combined <- ggplot(data = risk.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_continuous(limits = c(0,10), labels=c(0, 2.5, 5.0, 7.5, 10)) +
  geom_ribbon(aes(x=risk.work.dest, ymin=work.risk.potmig.lowerx, ymax = work.risk.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=work.risk.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=risk.work.dest, ymin=enter.risk.potmig.lowerx, ymax = enter.risk.potmig.upperx), color="grey", alpha=.2) + 
  geom_line(aes(x=risk.work.dest, y=enter.risk.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.risk.potmig.jitter),sides="t", color="grey") +
  geom_rug(aes(x=work.risk.potmig.jitter),sides="b", color="black") +
  xlab("Perceived Risk") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+ 
  theme(legend.title = element_blank())

risk.potmig.plot.combined

### Figure 4 ###

### Punishment and ener (illegal migration)
## running the model

enter.punishment.potmig <- lm(enterlist.pre ~ enterlist.pre.dum + punishment.enter + enterlist.pre.dum*punishment.enter
                              + gender + age + education + income.imp + freq.help.net + factor(parish), 
                              data= potmig)

summary(enter.punishment.potmig)

## generating clustered standard errors
enter.punishment.potmig.clust<-cluster.vcov(enter.punishment.potmig, potmig$ed.parish)
enter.punishment.potmig.mcse <- coeftest(enter.punishment.potmig, enter.punishment.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
enter.punishment.potmig.b2 <- enter.punishment.potmig.mcse[2,1]
enter.punishment.potmig.b3 <- enter.punishment.potmig.mcse[3,1]
enter.punishment.potmig.b22 <- enter.punishment.potmig.mcse[22,1] # interaction term
# variances
enter.punishment.potmig.varb2 <- enter.punishment.potmig.clust[2,2]
enter.punishment.potmig.varb3 <- enter.punishment.potmig.clust[3,3]
enter.punishment.potmig.varb22 <- enter.punishment.potmig.clust[22,22]
# covariances
enter.punishment.potmig.covb2b22 <- enter.punishment.potmig.clust[2,22]
enter.punishment.potmig.covb3b22 <- enter.punishment.potmig.clust[3,22]

list(enter.punishment.potmig.b2, enter.punishment.potmig.b3, enter.punishment.potmig.b22, enter.punishment.potmig.varb2, 
     enter.punishment.potmig.varb3, enter.punishment.potmig.varb22, enter.punishment.potmig.covb2b22, enter.punishment.potmig.covb3b22)

## transforming dataframe to graph    
enter.punishment.potmig.df <- potmig %>% mutate(punishment.enter = seq(1, 4, length.out=802),
                                                enter.punishment.potmig.conbx = enter.punishment.potmig.b2+enter.punishment.potmig.b22*punishment.enter,
                                                enter.punishment.potmig.consx = sqrt(enter.punishment.potmig.varb2 + enter.punishment.potmig.varb22*(punishment.enter^2)+
                                                                                       2*enter.punishment.potmig.covb2b22*punishment.enter),
                                                enter.punishment.potmig.ax = 1.96*enter.punishment.potmig.consx,
                                                enter.punishment.potmig.upperx = enter.punishment.potmig.conbx + enter.punishment.potmig.ax,
                                                enter.punishment.potmig.lowerx = enter.punishment.potmig.conbx - enter.punishment.potmig.ax)


### Punishment and work (semi-legal migration)
## running the model

work.punishment.potmig <- lm(worklist.pre ~ worklist.pre.dum + punishment.work + worklist.pre.dum*punishment.work
                             + gender + age + education + income.imp + freq.help.net + factor(parish), 
                             data= potmig)

summary(work.punishment.potmig)

## generating clustered standard errors
work.punishment.potmig.clust<-cluster.vcov(work.punishment.potmig, potmig$ed.parish)
work.punishment.potmig.mcse <- coeftest(work.punishment.potmig, work.punishment.potmig.clust)

## Pull out data for marginal effects plots
# coefficients
work.punishment.potmig.b2 <- work.punishment.potmig.mcse[2,1]
work.punishment.potmig.b3 <- work.punishment.potmig.mcse[3,1]
work.punishment.potmig.b22 <- work.punishment.potmig.mcse[22,1] # interaction term
# variances
work.punishment.potmig.varb2 <- work.punishment.potmig.clust[2,2]
work.punishment.potmig.varb3 <- work.punishment.potmig.clust[3,3]
work.punishment.potmig.varb22 <- work.punishment.potmig.clust[22,22]
# covariances
work.punishment.potmig.covb2b22 <- work.punishment.potmig.clust[2,22]
work.punishment.potmig.covb3b22 <- work.punishment.potmig.clust[3,22]

list(work.punishment.potmig.b2, work.punishment.potmig.b3, work.punishment.potmig.b22, work.punishment.potmig.varb2, 
     work.punishment.potmig.varb3, work.punishment.potmig.varb22, work.punishment.potmig.covb2b22, work.punishment.potmig.covb3b22)

## transforming dataframe to graph    
work.punishment.potmig.df <- potmig %>% mutate(punishment.work = seq(1, 4, length.out=802),
                                                work.punishment.potmig.conbx = work.punishment.potmig.b2+work.punishment.potmig.b22*punishment.work,
                                                work.punishment.potmig.consx = sqrt(work.punishment.potmig.varb2 + work.punishment.potmig.varb22*(punishment.work^2)+
                                                                                       2*work.punishment.potmig.covb2b22*punishment.work),
                                                work.punishment.potmig.ax = 1.96*work.punishment.potmig.consx,
                                                work.punishment.potmig.upperx = work.punishment.potmig.conbx + work.punishment.potmig.ax,
                                                work.punishment.potmig.lowerx = work.punishment.potmig.conbx - work.punishment.potmig.ax)


## interaction plot with confidence intervals based on clustered SEs

# The rug
enter.punishment.potmig.temp <- potmig$punishment.enter
work.punishment.potmig.temp <- potmig$punishment.work
a<-runif(802,min=0,max=0.6)
a<-a-0.2
enter.punishment.potmig.jitter<-enter.punishment.potmig.temp+a
work.punishment.potmig.jitter<-work.punishment.potmig.temp+a

# The plot
punishment.potmig.df <- cbind(work.punishment.potmig.df, enter.punishment.potmig.df)
colnames(punishment.potmig.df) <- make.unique(names(punishment.potmig.df))

punishment.potmig.plot.combined <- ggplot(data = punishment.potmig.df) + 
  theme(text = element_text(size=30)) +
  scale_y_continuous("Estimated Support") + 
  scale_x_discrete(limits = c(1,4)) +
  geom_ribbon(aes(x=punishment.work, ymin=work.punishment.potmig.lowerx, ymax = work.punishment.potmig.upperx), alpha=.2) + 
  geom_line(aes(x=punishment.work, y=work.punishment.potmig.conbx, linetype = "Semi-Legal"), size = 1) + 
  geom_ribbon(aes(x=punishment.enter.1, ymin=enter.punishment.potmig.lowerx, ymax = enter.punishment.potmig.upperx), alpha=.2) + 
  geom_line(aes(x=punishment.enter.1, y=enter.punishment.potmig.conbx, linetype = "Illegal"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=enter.punishment.potmig.jitter),sides="t", color="grey") +
  geom_rug(aes(x=work.punishment.potmig.jitter),sides="b", color="black") +
  xlab("Perceived Severity") +
  ggtitle("Aspiring Migrant Support\n for Unauthorized Strategies") + 
  theme(plot.title = element_text(face="bold", hjust = 0.5))+
  theme(legend.title = element_blank())

punishment.potmig.plot.combined

### Table 2 ###

### Model 1: enter (illegal migration) 

## running the model
enter.potmig <- lm(enterlist.pre ~ enterlist.pre.dum  + fair.jam + norm.family.rev
                   + enterlist.pre.dum*fair.jam +  enterlist.pre.dum*norm.family.rev + risk.enter.dest
                   + risk.enter.dest*enterlist.pre.dum + punishment.enter + punishment.enter*enterlist.pre.dum 
                   + gender + gender*enterlist.pre.dum + age + age*enterlist.pre.dum 
                   + education + education*enterlist.pre.dum + income.imp + income.imp*enterlist.pre.dum 
                   + freq.help.net + freq.help.net*enterlist.pre.dum + factor(parish), data= potmig)
summary(enter.potmig)

## generate clustered standard errors 

enter.potmig.clust<-cluster.vcov(enter.potmig, potmig$ed.parish)
enter.potmig.mcse <- coeftest(enter.potmig, enter.potmig.clust)
enter.potmig.mcse

### Model 2: work (semi-legal migration) 

## running the model
work.potmig <- lm(worklist.pre ~ worklist.pre.dum  + fair.jam + norm.family.rev
                  + worklist.pre.dum*fair.jam +  worklist.pre.dum*norm.family.rev + risk.work.dest
                  + risk.work.dest*worklist.pre.dum + punishment.work + punishment.work*worklist.pre.dum 
                  + gender + gender*worklist.pre.dum + age + age*worklist.pre.dum 
                  + education + education*worklist.pre.dum + income.imp + income.imp*worklist.pre.dum 
                  + freq.help.net + freq.help.net*worklist.pre.dum + factor(parish), data= potmig)
summary(work.potmig)

## generate clustered standard errors
work.potmig.clust<-cluster.vcov(work.potmig, potmig$ed.parish)
work.potmig.mcse <- coeftest(work.potmig, work.potmig.clust)
work.potmig.mcse
